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We study a homogeneously driven granular gas of inelastic hard particles with rough surfaces subject to 
Coulomb friction. The stationary state as well as the full dynamic evolution of the translational and rotational 
granular temperatures are investigated as a function of the three parameters of the friction model. Four levels 
of approximation to the (velocity-dependent) tangential restitution are introduced and used to calculate transla- 
tional and rotational temperatures in a mean field theory. When comparing these theoretical results to numerical 
simulations of a randomly driven mono-layer of particles subject to Coulomb friction, we find that already 
the simplest model leads to qualitative agreement, but only the full Coulomb friction model is able to repro- 
duce/predict the simulation results quantitatively for all magnitudes of friction. In addition, the theory predicts 
two relaxation times for the decay to the stationary state. One of them corresponds to the equilibration between 
the translational and rotational degrees of freedom. The other one, which is slower in most cases, is the inverse 
of the common relaxation rate of translational and rotational temperatures. 

I. INTRODUCTION 

Granular media are collections of macroscopic particles with arbitrary shape, rough surfaces, and dissipative interactions 
Many phenomenona are well reproduced by model granular media, where spheres are used instead of other, possibly 
more realistic shapes. In order to study such model systems, kinetic theories I4l l5l la. F7L l9l flOl fill fl2l fl3l fl4l ITM flal and 
numerical simulations 14 Ha. fl7l Ua. fl9l l20l l2lt l22l 12311 have been app lied for special boundary conditions and a variety of 
interesting experiments have been performed, see for example ll24l l25l l26l l27l l28ll . The dynamics of the system is usually 
assumed to be dominated by instantaneous two-particle collisions. These collisions are dissipative and frictional, and conserve 
linear and angular momentum while energy is not conserved. In the simplest model, one describes inelastic collisions by a 
normal restitution coefficient r only. However, surface roughness and friction are important fTol IT3I I20L I2TI I22I l29ll . since they 
allow for an exchange of translational and rotational energy and influence the overall dissipation. In the standard approach 
1I5 HI0I.I22I1 . surface roughness is accounted for by a constant tangential restitution coefficient r t , which is defined in analogy to r 
in the tangential direction. A more realistic friction law involves the Coulomb friction coefficient [i fl7ll30ll3lll32il . so that the 
tangential restitution r t (^) depends on the impact angle 7, i.e. the angle between the contact normal and the relative velocity of 
the contact points. 

Recently, Jenkins and Zhang lfl4ll proposed a kinetic theory for frictional, nearly elastic spheres in the limit of small friction 
coefficient /1. They introduced an effective coefficient of normal restitution by approximately relating the rotational temperature 
to the translational one. Thereby the kinetic theory for slightly frictional, nearly elastic spheres has the same structure as that 
for frictionless spheres. Also for small fi, Goldhirsch et al. 11611 showed that an infinite number of spin-dependent densities is 
needed to describe the dynamics of frictional spheres and that the distribution of rotational velocities is non-Gaussian. A mean 
field theory for three dimensional cooling systems of rough particles with Coulomb friction was proposed in fl3ll and found to 
be in very good agreement with computer simulations for a wide range of parameters. A systematic theoretical study of driven 
systems over the whole range of dissipation and friction parameters is not available to our knowledge. 

In the following, we propose a mean-field (MF) theory of homogeneously driven rough particles that accounts for Coulomb 
friction (i.e. a non-constant r t ) on different levels of refinement. The most accurate description parallels the three-dimensional 
(3D) results flUl for freely cooling systems. In addition, we present different levels of approximation to the full model and 
discuss their shortcomings in MF theory. The homogeneous driving used here is the same as in other recent studies of driven 
systems fllE^l . 

To test our analytical results we have performed numerical simulations of a randomly driven mono-layer of spheres, using an 
Event Driven (ED) algorithm l33ll . One key result is that, via ^(7), all parameters of the collision model affect the 

evolution of the translational and rotational degrees of freedom (temperatures) of the system. Only the full MF theory is able to 
quantitatively predict the system behavior for the whole parameter range. 
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The model system is introduced in section [H] The distribution of impact angles, as affected by translational and rotational 
degrees of freedom, is computed in section [HI] The standard approach with constant tangential restitution is briefly reviewed, 
before we introduce three levels of approximation and the full MF theory in section ITVl In section[V]we discuss the stationary 
state and and in section lVTl the dynamic evolution towards the stationary state. In both sections we compare the predictions of 
full MF theory and its approximations to simulations. Finally we present a summary and conclusions in section IVTT1 



II. MODEL 



The model system contains N three-dimensional spheres of diameter 2a, mass m, and moment of inertia I interacting via a 
hard-core potential. The particles are confined to a two-dimensional (2D) square with periodic boundary conditions. The linear 
box size is L and the area (volume) V = L 2 . The moment of inertia can be expressed using the shape factor 

o:=-^T- (1) 



For spheres with a homogeneous mass distribution q = 2/5. Inelasticity and roughness are described by a coefficient of normal 
restitution r, the Coulomb friction law with coefficient of friction \x, and a coefficient of tangential restitution r* which depends 
on r, fj,, and the impact angle 7 for sliding contacts, or on a maximum tangential restitution r™ for sticking contacts, when 
some "tangential elasticity" becomes important. In a collision of two particles i = 1 and 2 with positions contact normal 
n = (t*i — ?-2)/(2a), angular velocities a;, and relative translational velocity v\i = V\ — V2 (see Fig.^, their velocities 
after the collision are related to the velocities before the collision, through a collision matrix 12(1 13Z 13411 which is derived from 
the conservation laws for linear and angular momentum, energy/dissipation balance, and Coulomb's law. This three parameter 
model is able to reproduce the experimental measurements on colliding spheres of various materials 130113^1 . 




FIG. 1 : Schematic drawing of two-particle contact in the center of mass reference frame. Shown are the relative velocity g of the contact 
points, the impact angle 7 of the contact points, and the angle 712 between the relative translational velocity of the particles and their contact 
normal. 



A. Collision rules 



The collision rules are most transparent when written in terms of the relative velocity of the contact point in the center-of-mass 
reference frame 

g=t) 1 -D 2 -a(w 1 +w 2 )xn. (2) 

We decompose g = g n + g t into its normal and tangential components with respect to n, g n = (g ■ n)n and g t = g — g n . The 
change of normal momentum of particle 1, denoted by AP^™' is the same as for smooth particles 

ApW=-(m/2)(l + r) fln . (3) 
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The change of tangential momentum 

AP m = --^jm(l + r t )g t (4) 

is, in general, a function of the impact angle 7. Coulomb friction can be expressed l34ll in terms of a coefficient of tangential 
restitution 

r t (7)=min[rf( 7 ),rf] , (5) 

which is a function of the impact angle 7 between g and n. Here rj™ is the coefficient of maximum tangential restitution, with 
— 1 < < 1 to ensure that energy is not created. The quantity rf (j) is determined using Coulomb's law 

r? (7) = -1- 2-^-^(1 +r) cot 7, (6) 

q 

with the impact angle tt/2 < 7 < tt so that cos 7 = g ■ n/\g\ is always negative I20ll30ll32il . Here, we have simplified the 
tangential contacts in the sense that exclusively either Coulomb friction applies, i.e. AF' 1 ' = fiAP' n \ or constant tangential 
restitution with the maximum tangential restitution coefficient r™. Coulomb friction is effective when the relative tangential 
velocity is large, whereas tangential restitution applies for low tangential velocities. 

Note that in the general case v Iot = —a(u>i + u)%) x n / 0, so that the angle 712 between the contact normal n and the 
relative translational velocity v±2 = v\ — v-x is different from the impact angle 7 of the contact points, see Figs.QJand|21 In the 
following we will refer to 7 when we talk about the impact angle. 



1.0 




FIG. 2: Tangential restitution rt as function of the impact angle 7 for different values of the coefficient of friction fj,. 



B. Driving model 

The drivin g of a g ranular material can be realized by moving walls, see Ref. GJ and references therein, corresponding to a 
local heati ng [3^[37ll38ll. or the s yste m can alternatively be driven by a global homogeneous, random energy source in different 
variations lllll |12 

HUMES Eli 

We choose homogeneous translational driving here and modify the velocity of particle i 

at each time of agitation t such that 

v'^t) = Vl (t) + v dl &{t) (7) 

where the prime on the left hand side indicates the value after the driving event. Measuring masses in units of the particle mass 
to, the driving velocity v^t sets the time (velocity) scale and defines the driving temperature Tdr := muL The components of 
the vector £j(t), £i, x (t) and fj^Ct), are uncorrected Gaussian random numbers with zero mean and variance 



(6,k(*)&i(0>tt} = MM*(t-0. 



(8) 
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where 5{j and 5ki are Kronecker deltas and S(t — t') is the Dirac delta function. The stochastic driving rule in Eq. leads to an 
average rate of change of temperature 

AT/ At = H dl , with H dl = f dl T dr , (9) 

after every driving time-step At = f^ T . 



C. Simulations 



We have performed simulations of a randomly driven mono-layer of sphere s, using an Event Driven (ED) algorithm 12(1 EH 
I29H43I1 . and compared the results with the MF predictions, see also Refs. ITU [13I.I29T l4ll l42ll . Every simulation is equilibrated 
without driving with r = 1 and in the smooth surface limit r" 1 = — 1 . Then inelasticity, friction and drivi ng a re switched on, 
according to the rules defined above. The problem of the inelastic collapse characteristic of the ED alg orithm 14411451 . is handled 
by using normal restitution coefficients dependent on the time elapsed since the last event 13 ■ The frequency of driving 

is chosen such that it is larger than or comparable to the typical collision frequency per particle, both initially and in steady state. 
Varying the driving frequency to much larger values did not affect the simulation results, whereas the use of a much smaller 
driving rate caused different results due to the slow input of energy. 



III. IMPACT-ANGLE PROBABILITY DISTRIBUTION 



In the following we shall discuss various levels of approximation to the collision rules given in Eqs. © and (|6jl. One possibility 
to simplify the collision rules is to consider tangential restitution averaged over all impact angles 7, thereby reducing the problem 
to one with a constant coefficient of tangential restitution. For that purpose we need to know the probability distribution of impact 
angles. 

The assumption of "molecular chaos" implies a homogeneous distribution of the collision parameter b = 2a sin 712 which 
is simply related to the angle 712 between the relative translational velocity V\2 and the contact normal n according to 
cos 712 = V12 • / 1 1» 12 1 , see Fig.^ Hence the probability distribution of sin 712 is constant, P{ 2 (sin.7i2) = 1. (The "prime" 
indicates probability functions of the sine or the cosine of the angle.) A uniform probability P' implies for the distribution of 
the angle ^12(712) = — cos7 12 , so that grazing contacts appear less probable than central collisions when a fixed interval d7 12 
is considered. The uniform P[ 2 (sin 712) is in agreement with our numerical data, see Fig. [3] 




FIG. 3: Plots of the probability distribution of 7 from simulations (symbols) and from Eq. 1121 with R values from the simulations. The arrows 
indicate the corresponding 70, while the parameters are (a) r — 0.95, n = 0.5 and variable r™, and (b) r = 0.95, r™ = 0.4 and variable ^t. 

In general, the impact angle 7 between the relative velocity of the contact point g and the contact normal n is different from 
the angle 712 between the relative translational velocity V12 and the contact normal n, as displayed in Fig. [2 The two angles 
are identical only in the case of smooth particles or in the limit of vanishingly small rotational velocities. In the general case we 
compute P'(cos7) by averaging over all binary collisions 

^H'H-isf)).,,- <10) 



5 



This average can only be computed approximately. We assume that the translational and rotational velocities of the colliding 
particles are distributed according to Gaussians with a temperature T tr for the translational and a temperature T rot for the 
rotational velocities. Within this approximation the above average is given explicitly by 



P'(cos 7 ) 



J IS 



cos 7 



gn 



J(l) 



(11) 



with the phase space integral 



J{X) = J dY 1 dT 2 («ia ■ n) 0(-«! 



n) S(\r 12 \-2a)X , 



where X = X(Ti,T3), and the phase space element 

dT k = d\ k d 2 v k <Lj k e-™i/( 2T ^e- Iw i/V T ^ 

for k = 1,2. 

The remaining integrals can be computed analytically, yielding the following expression for the impact angle distribution 



P{l) = - 



{1 + R/q) cos 7 
(1 + [R/q] cos 2 7 



,3/2 



(12) 



Here we have introduced the ratio of rotational and translational temperatures R := T Iot /T tI and recall q = I /(ma 2 ). The 
probability distribution P(j) is compared to the results of our simulations in Fig. [3] reasonably good agreement is observed. 
With increasing rotational velocities, contacts with large g t (small 7) become more and more frequent due to the increasing 
rotational contribution. On the other hand, collisions with vanishing g t (large 7) become less probable, since the rotational 
contribution leads to a net increase of g t . 



IV. DIFFERENTIAL EQUATIONS IN MEAN FIELD THEORY APPROXIMATIONS 

In the following we present different approximations for frictional particles, referred to as models A-E. Model A is the well 
known model using constant coefficients of normal and tangential restitution, cf., e.g., fl fioll . Model E implements Coulomb 
friction as introduced by Walton fl7ll . While model A is the mean field solution for rough particles with a constant coefficient 
of tangential restitution, model E is the mean field solution for particles with Coulomb friction. Models B through D are 
approximations to model E that may be simpler to deal with but have significant shortcomings. 

The starting point of our mean-field approach is the theory of Ref. Holl for a freely cooling gas of rough particles with a 
constant coefficient of tangential restitution (r t = const., corresponding to the limit fx — > 00). The theory is based on a 
pseudo-Liouville-operator formalism and on the assumption of (i) a homogeneous state, (ii) independent Gaussian probability 
distributions of all degrees of freedom, i.e. all components of the translational and the rotational velocities, and (iii) the assump- 
tion of "molecular chaos", i.e. subsequent collisions are uncorrected. The agreement with simulations is very good as long as 
the above assumptions are valid I2H1 . 

The main outcome of this approach is a set of coupled time evolution equations for the translational and rotational MF 
temperatures T tr and T Iot fioll which can be extended to also describe arbitrary energy input (driving) fl5l l29l E2I1 . Given the 
random driving temperature Tdr and an energy input rate /d r , as defined above, one just has to add the positive rate of change of 
translational energy H^r, see Eq. (|9}, to the system of equations 12911 . 



A. Model A: Constant tangential restitution r t = r™ 

We recall the results of the mean field theory for the model with a constant coefficient of tangential restitution which is 
obtained from the general case in the limit \i — > 00 (see Eqs. (14) in Ref. 12 110 . The system of coupled equations reads in 2D: 

5t^tr(*) = #dr + G 

£T rot (t) = 2G 

Note the choice of signs which lead to positive coefficients. Based on more physical arguments, A quantifies the dissipation 
of translational energy, B and B' correspond to the interchange of energy between the translational and rotational degrees of 



-ATI' 2 



B'T t 



3/2 



f BT t \ /2 T ro 

^ - 1 tr 1 rot 



(13) 
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freedom, and C describes the dissipation of rotational energy. The coefficient G sets the time-scale of the system, i.e. the 
collision rate (per particle) r _1 = (l/2)GT t 1 / 2 , with 

G=-^=vg 2a {p). (14) 

Here g2aiy) denotes the pair correlation function at contact. In the approximation proposed by Henderson ll5l l49ll50H5lll52il . 

Qiaiy) = (1 — 7z//16)/(l — v) 2 , it depends only on the 2D volume fraction of the granular gas v = na 2 N/V. The four constants 
A, B, B' and C read in this limit 

1 -r 2 

A = A r + Ajj a , A r :— — - — , (15) 
A V0 := f(l-%), (16) 

2 

B' = B = B V0 := and (17) 



It is useful to define a function 



^ : =^T^' fOT0 ^^^ <1 - (19) 



which has to be evaluated at constant tangential restitution r t = r" 1 in the limit /x — ► oo 



B. Model B: Simplified mean tangential restitution n — (rt) 12 

A first step beyond the above theory with a constant r/o = V^t 1 )^ i s tne replacement of 7^(7) by its average 

(rt) - J d 7 P(7)rt( 7 ). (21) 

7T/2 

The integral over 7 from 7r/2 to 7r, has to be split into two parts, one corresponding to the range it/ 2 < 7 < 70 for which there is 
Coulomb sliding with r t given by Eq. 0, and a second part corresponding to the range 70 < 7 < it, for which there is sticking 
with constant r t = r" 1 (see Fig.|2)- The critical angle 70 is given by 

all + r m ) 

c:=-cot 7o = ( q ^*' ; >0- (22) 
To simplify the computation, we use the approximation P(7) ~ ^12(7) = — cos (7), such that 



9+1 

\' t/12 — — n 

with the abbreviation 



r t )i 2 =-l + 2±i(l + r)/i ln(c + /) . (23) 



/ := y/l + c 2 . (24) 

The averaged coefficient of tangential restitution (r t )i2 must be inserted into 77 in Eq. dl9l l. Thus we obtain the same set of 
coefficients as in Eqs. dl5l-dl 81 with ?/ replaced by 

Vi :=v((rt) 12) = -ln(c + f) . (25) 

c 

In this approach, only the average value of r t is considered and fluctuations of ?' t with 7 are neglected. Furthermore the 
difference between 7 and 7x2 has been ignored in the averaging procedure. In contrast to model A this is the simplest model to 
incorporate the coefficient of Coulomb friction [i, (r t ) 12 = (?'t)i2(A t )- 
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C. Model C: Mean tangential restitution r t = (rt) (R) 



In model C we again replace rt (7) by its average but use the correct impact angle probability distribution function P(j) from 
Eq. i ll 2t in the averaging procedure. The result is an independent averaged coefficient of tangential restitution 



/ 1(nl -, q + 1 1 + r fi 

j [ f(/-c) 2 (x/-/ + cf) 
\(x/-/-cf)2(x/+/-cf) 



(26) 



with 



= x\R) := 1 + R/q, 



(27) 



/ = /(i?) := Vl + zV, (28) 

and / defined in Eq. i24\ . Note that x is an implicit function of time through R. For R — > (x — > 1) Eq. J26l > reduces to Eq. 
(1231 - as expected. For i? — > 00 (a; — > 00) there is no friction and (rt)(R) — > — 1. 

We formally get the same differential equations dl 31 but with non-constant coefficients A = A(R), B' = B = B(R), and 
C = C(R) which are obtained by replacing 770 by r)((r t ) (R)) in Eqs. J15b - dl8l >. These coefficients are implicitly time dependent 
via R. 



1. Constant tangential restitution limit 



In the limit fi — > 00, c — > 0. In that case model C reduces to model A. 



2. Weak friction limit 



For /x — > 0, c — * 00 we recover smooth spheres with (rt) — >• — 1. A series expansion to lowest order in (equivalent to lowest 
order in c" 1 ) of Eq. i26\ reads 

(r* ) (R) = -1 + — (1 + r) £ { I In ( M ) | + In (x)+ 
q X I 

ln (r^)} +G(Al3) ' (29) 

expressed in terms of x and /i. 

As long as x stays finite (which is the case for a driven system) the leading order is thus fj,\ In (/i) | for small /x. For x — > 1, Eq. 
( I29> yields the same result as Eq. d23i in leading order in c _1 . 



i. Comparison of model B and model C 



Due to the implicit nature of model C it is rather difficult to work out its predictions, e.g., for the ratio of temperatures. 
Therefore, we present here the mean tangential restitution from models A, B, and C in Fig.|4] Note that (rt) for model C 
depends not only explicitly on [1 but also implicitly through R. To keep the discussion simple, we present results only for some 
constant, representative values of R. The mean restitution for large R is smaller (or equivalently, the corresponding /j, is larger) 
than for small R. Models B and C become indistinguishable in the limit R — > 0, as expected. 
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FIG. 4: Expected mean tangential restitution, (rt), as function of the friction coefficient fi for models A, B, and C. The parameters used are 
r = 0.95, rT = 1.0 (for A, B) and different R = 1.0, 0.40, and 0.15 (model C: solid lines from right to left). 

D. Model D: Variable (simplified) tangential restitution r t (712 ) 

In this section and in the following one, we discuss a coefficient of tangential restitution which depends on 7. Model D is 
defined by approximating 7 sw j r2 , which is strictly true only for R — ► or fj, — > (or equivalently r™ — ► —1). We again obtain 
the same differential equations Jl 31 for T TO t and Tt r with the coefficients 



A = A, = A r + [A vo + A*} /f 
B = = [B m 



(30) 



B' = B> = [B 



no 



B*]/f, 
B'*} If . 



C=C^ = [C V0 + C*\ If , 

and / defined in Eq. J24l >. The terms that originate from Coulomb sliding are denoted by an asterisk and are given explicitly by 



B* = 

B" = 
A* = 
C* = 



9 9 



2q(f + lf ' 
(2f + l)B* , 
7yoc 2 /2 — qB* , and 
( Vl f - ?/0 - 2B*) l(2q) , 



(31) 



expressed in terms of / [cf. Eq. (I24H . 770 [cf. Eq. J20H . 771 [cf. Eq. (I25H . and q [cf. Eq. Q]. The terms B* and B'* are strictly 
positive, while the dissipation correction terms A* and C* , in principle, can change sign. Note also that B* and B'* are not 
identical here. All coefficients depend on the system parameters only. They are constants in time - in contrast to model C (and 
E as will be shown later). 



1. Constant tangential restitution limit 



In the limit \i — ► 00, one has c — > 0, i.e. / — > 1, and all correction terms {A* , B* , B'* , C*\ — > so that one obtains Eqs. 
dl 31 -fll 81. Note in particular that the coefficients B^ and B'^ are equal only in the limit /i — > 00 iflOtl . 
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2. Weak friction limit 

In the limit /.t — * (c — * oo, / — * c), the lowest order expansion in c _1 leads to an approximation of the coefficients in Eqs. 
( I30> . where we have used i] /c = fx(l + r) /2: 

B, = -^V + Oiv 2 ), (32) 



q 2 



1 fl + r N 2 



3. = ll^T 1 ) ^ 2 +0(M 3A 



9 V 2 

1 + r 



= A r + — - — fj, + 0(fi 

,-±l±r,(| bM i +ta (^L)-,2 

+ 0(e 2 ) ■ 

From Eqs. d32i . we learn that B'^ is second order in whereas B M is first order in /i, reflecting an asymmetry in the energy 
transfer rates. On the other hand, A^ as A r is almost constant, whereas C M depends on \i logarithmically which is an artifact of 
our approximation 712 ~ 7, see Eq. d35l below. 



E. Model E: Variable (exact) tangential restitution ^(7) 



The final step of refinement of the MF theory is to use 77(7), instead of 7 , t (7 12 ), to compute the coefficients. This is the full 
mean field theory. The calculation is similar to the one for 3D in JT3I1 and is presented in appendix|X] We obtain the following 
coefficients, to be inserted into Eqs. (1131 . 



A = A fl (R) = A r 



A 



B = B„(R) = 
B'=%{R) = 
C = C^R) = 



B 



B,, 



no 
B* 

B' 



C?)o "F C 



A* 
I? 

A 3 
// 3 



7 3 



(33) 



(34) 



with /, x and c defined in Eqs. (I28i . d27t and (I22> . respectively. The new correction terms are in detail: 

B* = -Voc 2 /(2q), 

(2/ + l)(ry cx 2 ) 2 
2g(/ + l)2 

A* = -q (^B* + B'*^j , and 
C* = -x 2 B* , 

with q and x as introduced in Eqs. Q and i21\ . Interestingly, we find now a negative B* together with positive coefficients B'* 
and C*; only A* can be both positive and negative. Like in model C but in contrast to models A, B, and D, here the coefficients 
are implicit functions of time, again. 

In conclusion, models D and E appear similar in shape but there are several striking differences: (i) The division by / and 
/ 3 in model D is in contrast with the division by / 3 in model E, (ii) the term B* in model D is always positive, while B* in 
model E is always negative, (iii)the sign of C* in model D is not determined a-priori, while the term C* is always positive, (iv) 
among the correction terms of model E, only B* is independent of R, and (v) the more refined theory appears in a simpler form, 
especially the term C*. 



1. Constant tangential restitution limit 

The limit of constant tangential restitution can be reached by taking the limit /j, — > 00. In this case c — > 0, / — ► 1 and thus all 
additional coefficients A*, B*, B'*, and C* vanish such that Eqs. J 1 3I >- J 1 8I > are recovered. 
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2. Weak friction limit 

In the limit \i — > (c — * oo, / — + xc) an expansion to the lowest order in fi leads to an approximation of the coefficients in 
Eqs. (1331 when we remember that r/ /c = (1 + r)^i/2: 

MR) = -^^ + (^ 3 )' < 35 > 



A^R) = A r -qf K B fl (R)+B^R))+0{ f i i 
C„{R) = ~x 2 B^R) + 0(^) . 



Since x = x(R) approaches one in the weak friction limit, both B^(R) and C fJi (R) are proportional to [i in leading order. To 
lowest order in fi, Eq. (1351 predicts A^(R) = A r + O(fi), i.e. proportional to while B' (R) is proportional to /z 2 . 
For Eqs. (II 31 with J35i simplify to 

jT t r{t) = H dT - GT t 3 / 2 (i^ + 0(jij) , (36) 

which means that in the limit of low friction the differential e quat ions for T tr and T rot decouple. In the non-driven case this 
leads to surviving rotational energy (not show), similar to Refs 



quations 



V. STEADY STATE 

Before discussing the approach to the stationary state in the next chapter, we first elucidate the stationary state and compare 
results of our simulations to various levels of refinement of the mean field theory. 

A. Analytical results 

By imposing ^T t s * at = and ^T r s * at = one gets the steady state values of the rotational and the translational temperatures. 
For models A, B and D, the coefficients in the differential equation do not depend on R (or x). Therefore the solution is simply 

T(TT \ 2/3 
~gi ) ' 

with 

R stat = B'/C , and X = A — BR stat , (38) 
as discussed in more detail for all models in the following. 

/. Model A 

For model A, the steady state ratio of rotational to translational energies is 



and the energy dissipation factor is 



E stat = (39) 

q-vo 



4 2 2(g-77o) 



Note here again that model A does not contain any dependence on the coefficient of friction \i. 
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2, Model B 



Model B evolves from model A, by just replacing 770 by 771 (/j) = (770 /c) In (c + /) from Eq. d25i in the above two Eqs. ( I39t 
and l !40t . so that, e.g., 

^stat = g ? 7i = gfao/c) In jc + f) 

q-T)\ q- (770/c) In (c + /) ' 

In the limit of small /jCI, the leading order terms are JJ stat « (1 + r)/i| In /i|/2 and X « i^ 1 — h In 



3. Morfe/ £> 



From model D, the following, more complex terms are obtained: 
and 

so that, asymptotically for /j< 1, model D leads to behavior similar to that of model B. 



4. Model E 



Formally, we can write down Eqs. fl37i for model E, too. Instead of using Eqs. j38t . i? stat must be extracted (numerically) 
from Eq. JA22> where the left hand side vanishes in the stationary case. It can be show analytically that there is always a unique 
solution - in contrast to the freely cooling case fl3ll . With the solution for R stat at hand, Eq. dA19l > (with a vanishing left hand 

side) can be written in the form T t s r tat = (%y) 2 ^ a S a i n where I is a nonlinear function of i? stat whose particular form can be 
easily seen from Eq. ( IA19I . 



5. Models C and Efor small fi 

For models C and E the coefficients in the differential equations do depend on R, so that the steady state values have to be 
computed numerically for a general choice of parameters. Analytical results can only be achieved in the limit \i -C 1, where we 
can use the expansions of the coefficients introduced in sections llV Cl and flV El 

For model C we obtain to lowest order in /i, the dissipation factor T « A r and, using 772 = (??o/2) In c, 

m«i cy \cqj q\cJ m<i 2 
For model E we find again I w A r and 




very similar in shape to the result from model C, besides the logarithm lnc that is hidden in the definition of 7/2- This leads 
to the qualitative difference in asymptotic behavior between models C and E: The correct asymptotic behavior for small fi is 
i? 8tat oc /_t. Note again that the more refined model E leads to a simpler analytical result than the approximated model C. 
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6. Discussion 

The expansions for small fj, <C 1 show that the result for i? stat based on model E, see Eq. J42l >. disagrees with all other models. 
In model E we find that i? stat vanishes linearly as ri — > 0, whereas models A-D predict a slower decrease, encoded in the \i\ In p\ 
dependence. Models A and B have the same analytical form for i? stat if expressed in terms of ?/ for model A and in terms of rji 
for model B. Similarly, models C and E have the same functional dependence on rj, if 772 is used for model C and 770 f° r model 
E. The comparison of the models for arbitrary values of p will be given in the next subsection, where we also present the results 
of our simulations and compare them to the predictions of the various mean field models. 

B. Comparison with simulations 

In this subsection, the steady state predictions from our models are confronted with the numerical simulation results. Note 
that we present results for rather high densities and dissipation, where our assumptions about homogeneity of the system and the 
Gaussian shape of the velocity distributions is not strictly true anymore. However, we want to stress the point that the present 
theory is astonishingly close to the numerical simulation with experimentally relevant parameters even when the most basic 
assumptions are somewhat questionable. 

1. Variation ofr™ 

In Figs.|5](a-c), the stationary rotational and translational temperatures and their ratio R are compared for r = 0.95, [i = 0.5 
and different values of r" 1 ; note that the data in (a) and (b) are scaled with the expression for /1 = 0. The symbols correspond 
to simulation data, with the error bars showing the standard deviation from the mean values. The lines correspond to different 
refinements of the theoretical approaches, i.e. models A, B, D, and E. 

For r' t n ~ — 1, the simulations agree with all theoretical predictions; for r" 1 w 1, large discrepancies are evident. The more 
refined a model used, the better the quality of agreement. The qualitative behavior of the data is best captured by model E, and 
we relate the remaining quantitative deviations to the fact that the simulations involve rather high density v and comparatively 
strong dissipation r. 

2. Variation of /i - translational temperature 

In Fig.[6]we plot the translational temperature in the same way as in Fig.[5Ja), but now, we keep the values r™ = 0.4 (a) and 
r™ = 1.0 (b) fixed and vary it. Furthermore, we compare data for r = 0.99 and r = 0.95 in one plot and observe satisfactory 
agreement between simulation results and the full mean field theory, model E. (The predictions from models A and B are only 
shown for r = 0.99.) 

For (realistic) values of r™ = 0.4, see Fig.|6ja), one obtains a transition from the \i = limit to the p, — > 00 value of the 
kinetic energy, over three orders of magnitude in p, whereas for r r " = 1.0, see Fig.|6jb), the kinetic energy first decays with p 
but then increases again to the stationary state temperature of smooth particles, since no energy is dissipated due to tangential 
friction for p — ► 00 and r™ = 1.0. 

Here, we remark that model A, with r t = r" 1 and the limit p — ► 00 is inadequate to model the ^-dependency of the data, 
it only gives the p — > 00 limit, as expected. Approach B only shows qualitative agreement with our simulation data, whereas 
theory D shows good quantitative agreement for small p. The agreement seems better for weak normal dissipation r = 0.99, as 
compared to the cases with r = 0.95. The deviations between simulations and model D in the intermediate range of p are due 
to values of R of the order of unity, for which the assumption 712 ~ 7 is not true, as pointed out above. 

For weaker normal dissipation r, one obtains a stronger reduction of the translational temperature in the range of strongest 
total dissipation (around p w 0.4). This is due to the comparatively stronger contribution of tangential dissipation. However, as 
in the previous subsection, the agreement between simulations and model E is satisfactory, especially for r — ► 1. 

3. Variation of p - rotational temperature 

In Fig.0we plot the ratio of rotational and translational temperature in the same way as in Fig. [3c), but now, like in Fig. [6] 
we keep the values r™ 1 = 0.4 (a) and r™ = 1.0 (b) fixed and vary p. Also here, we compare data for r = 0.99 and r = 0.95 in 
one plot. For the values of r™ examined (see Fig.0 one observes a smooth transition of R over about three orders of magnitude 
in p., from the value R = (in the limit p = 0) to the value R = r" 1 (in the limit p — » 00). Note that the observation R = rj" 
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FIG. 5: Simulation results (symbols) and theory (lines) for the parameters v — 0.34, N = 11025, r = 0.95, and ^ = 0.5, plotted against the 
maximum tangential restitution r™. (a) Translational temperature T t s r tat , and (b) rotational temperature T^t , plotted against r™, and scaled 
by T t s r tat (/i = 0), the mean field value for smooth particles, (c) Ratio of rotational and translational temperature 7?, plotted against r™. 



is coincidence, since the correct asymptotic result for large \i is R = 2(1 + r r t n )/(9 — 5r™). Again, the agreement between 
simulations and model E is impressive. 

All models agree qualitatively in the large /i-limit, even though the quantitative agreement with simulations is again best 
caught by model E, as can be seen in Fig. [8] 

The remaining question is the asymptotic behavior for very small fj,, as can be viewed in Fig.|5] and as discussed theoretically 
in subsection lV Al The quantitative behavior of R for small fi is tested by a power law fit of the numerical values, according to 
an expression R = b^ a . The fit gives a = 1.00(4), for r = 0.99, rf = 0.4, 1.0 and a = 0.99(4), for r = 0.95, r t m = 0.4, 1.0. 
Thus the asymptotic behavior is proportional to /i, in excellent qualitative and quantitative agreement with the prediction of 
model E. 
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FIG. 6: Translational temperature T t s r tat scaled by the mean field value for smooth particles T t s r tat (fi = 0), plotted against fi, for the parameters 
as in Fig.[5] The tangential restitution coefficients are fixed to (a) r™ = 0.4, and (b) r™ = 1.0. Data with normal restitution r = 0.99 (solid 
symbols and thick lines) and r = 0.95 (open symbols and thin lines) are compared. 

VI. APPROACH TO STEADY STATE 
A. Close to steady state 

Provided the system is sufficiently close to steady state, we can linearize the set of Eqs. Jl 3i around T t s r tat and T, s * at . This is 
particularly simple for models A, B, and D, where the coefficients in the differential equation do not depend on R and hence can 
be solved analytically for the stationary state. We set T tI (t) = T t s r tat (l + 5T tl (t)) and T rot (t) = T r s * t at (l + ST mt (t)) and obtain 
the linearized dynamic equations 

^r u . = G rr{(|4 + ^>r,, + ^r4 

^-6T Iot = 2GCT^ {ST tI - ST Iot } . (43) 
at 

This set of linear equations is easily solved to yield two relaxation rates Ai and A2. In a stable stationary state they must be 
positive and they are. We present here only results for the simplest model (A) and postpone the general discussion to the next 
paragraph, where the full dynamic evolution towards steady state will be examined. 

In Fig. [TO] we plot the two relaxation rates as a function of r™ for a fixed value of r = 0.95. In the limit of smooth spheres 
one of the rates vanishes because the rotational energy is conserved in that limit. For r™ ~ —0.84 the two rates are equal and 
for increasing r™ the difference between the two rates increases monotonically with r"\ such that for perfectly rough spheres 
the larger rate is about fourteen times the smaller one. Such a pronounced separation of time scales is familiar from the cooling 
dynamics of the same model, see lll3ll . There it was shown that the ratio of translational to rotational energy, R, relaxes fast to its 
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FIG. 7: R plotted against /Li, for the same parameters as in Fig.[6] The tangential restitution coefficients are again fixed to (a) j-J" = 0.4, and 
(b) r-r = 1.0. 
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FIG. 8: Deviation from equipartition, 1 — R, plotted against the inverse friction coefficient, fi 1 , for simulations from Fig. [fjb). Note the 
double-logarithmic scale of this plot. 

stationary value, whereas both the translational as well as the rotational energy decay on the same, much longer time scale. This 
point will be discussed in a more general setting (model E and relaxation from an arbitrary initial condition) in the subsequent 
paragraph. 

B. Full Dynamic Evolution 

In Fig. ^2 the full dynamic evolution of the translational and rotational temperatures with time is shown for two simulations 
with N = 11025, v = 0.0866, r — 0.95, r r " = 1.0, and different values for the coefficient of friction. In both situations, the 
agreement between simulations and the numerical solution for the full MF theory, model E, is good - not only concerning the 
limiting values and the asymptotes, but also the time dependence during the two regimes (i) equilibration between T tT and T rot , 
and (ii) approach to final steady state. 
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FIG. 9: Ratio of rotational and translational temperature, R, plotted against fi, for some simulations from Fig.|6fb). Note the double-logarithmic 
scale of this plot. 
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FIG. 10: Relaxation rates Xi : 2, close to steady state for r = 0.95 as a function of r™. 

VII. SUMMARY AND DISCUSSION 

In summary, a dynamic MF theory for the full time evolution of the translational and rotational temperatures of a homoge- 
neously driven two-dimensional granular gas has been presented. Particle collisions were modeled using the Walton model 11711 . 
i.e. with normal dissipation, tangential restitution (sticking) and Coulomb friction (sliding). The Walton model can be formulated 
in terms of a coefficient of tangential restitution, which depends on the impact angle 7. Using a Pseudo-Liouville operator we 
have computed the distribution of impact angles as well as the mean field dynamics and steady state values of the translational 
and rotational temperatures. 

In addition to the complete mean field theory of the Walton model ("model E"), we discussed three levels of approximation 
in order to simplify the differential equations of the time evolution. The crudest approximations including Coulomb friction 
("model B" and "model C") assume that an effective constant tangential restitution exists and can be computed by averaging 
over the angular distribution of impact angles. For model C this averaged coefficient depends on the current values of the 
translational and rotational temperatures and thus on time. Even simpler is model B where the rotational contribution to the 
impact angle is neglected, leading to a coefficient of tangential restitution that only depends on global system parameters. The 
closest approximation ("model D") to the full mean field theory ("model E") keeps the dependence of r t (7) on the impact angle 
7 but, like for model B, the contribution of the rotation of the particles to the impact angle is neglected. 

The predictions of the increasingly refined models of frictional dissipation as well as the full MF theory have been compared 
to simulations of a randomly driven mono-layer of spheres using an Event Driven algorithm. Emphasis has been put on the 
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FIG. 11: Evolution of temperatures with rescaled time, with r _1 = (l/2)GT tI (0) 1/2 , for simulations with TV = 11025, v = 0.0866, 
r = 0.95, r? = 1.0, and (a) fi = oo, (b) fi = 0.5. 

stationary state which is characterized by two temperatures, T tr and T rot , one for the translational and one for the rotational 
degrees of freedom. Guided by the MF approach we discovered a rich phenomenology like a non-trivial dependence of the 
stationary state temperatures on the model parameters. For example, the translational temperature is non-monotonic as a function 
of maximal tangential restitution r" 1 and also non-monotonic as a function of Coulomb friction jtt, provided r" 1 is sufficiently 
large. 

All models predict steady state values of the translational and rotational temperatures, which are considerably improved as 
compared to the model without friction ("model A"), which assumes constant tangential restitution (see Figs. [6] and Q. All 
approximations A-E agree in the limit of large friction, where the tangential restitution becomes independent of the impact angle 
(see Fig. [3. Qualitative agreement between models B-D and simulations is achieved also for intermediate values of /1. However 
in the limit fj, — ► all approximations break down and only the complete mean field solution ("model E") is in agreement with 
the simulations (see Fig.|9jl. In particular model E predicts the linear dependence of the ratio of temperatures, R = T Iot /T tr , on 
the friction coefficient fi that is observed in the simulations and was used in Ref. fl4ll to derive an approximate kinetic theory of 
frictional particles. 

Sticking contacts become more important relative to sliding contacts for fixed fi and decreasing r" 1 . In this regime models 
B and D seem reasonable, but lead to poor quantitative agreement as r"' approaches 1. The full mean field theory ("model E") 
leads to reasonable agreement for all values of r™. For weak dissipation, r — > 1, the agreement is very good - for stronger 
dissipation, we relate the deviations to the failure of both the homogeneity assumption and the molecular chaos assumption 
made. 

Linearizing the dynamic MF equations around the steady state leads to an eigenvalue problem with two relaxation rates, one 
of them being related to the equilibration between the translational and the rotational degrees of freedom, while the other one 
controls the approach of the system to its steady state. For strong coupling, the former process is much faster, so that there is a 
clear separation of time scales, which has been discussed already for a freely cooling system in the absence of driving. 
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In conclusion, realistic Coulomb friction turned out to be a subtle problem as only the full mean field theory of the Walton 
model predicts the effects of friction for all values of fi and r™. All simplifications are both qualitatively and quantitatively wrong 
in some parameter range. Our studies can easily be extended to three dimensional systems or more complex ones, like e.g. a 
polydisperse mixture of frictional particles with different material properties. Other driving mechanisms could be employed as 
well. 
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APPENDIX A: DERIVATION OF THE DIFFERENTIAL EQUATIONS FOR MODEL E 

The details of the derivation of the coefficients in Eqs. J33i and d34t for model E in Sec. llV El will be shown. The calculations 
are performed using a Pseudo-Liouville operator formalism fl3ll2lL l53ll . They are very similar to the ones in three dimensions 



11311 . First, we briefly recall the Pseudo-Liouville operator formalism. 

Let the vectors of position, translational and rotational velocity of a particle k, in a two-dimensional plane (x,y) with only 
vertical spin (z), be defined as r k = {r k , x ,r k . y ,0), v k = (v k , x , v k>y , 0), and u k = (0, 0, uj k ) . 

The time evolution of a dynamic variable A(t) that depends on time only through the positions and velocities of N particles, 
can be determined by means of a pseudo-Liouville operator C+ for t > 

A(t) = exp(iC+t)A(0) . (Al) 

The pseudo-Liouville operator £ + consists of three parts £ + = Co + C + + C*r. The last part, C*f, describes the homogeneous 
driving, the first one, Co, describes the free streaming of particles 

iY 



C = -i^2v k ■ V rfc , (A2) 



fe=i 



and the second one, C + = | Y^ijtk ^+ describes hard-core collisions of two particles 

Tf = i(v kl ■ r kl )G(-v kl ■ r kl )5(\r kl \ - 2a)[bf - 1). (A3) 

The operator bf replaces the linear and angular momenta of the two particles k and I before collision by the corresponding 
ones after collision, according to Eqs. and @. Q(x) is the Heaviside step-function, and we have introduced the notation 
r k i = r k — ri and r k i = r k i/\r k i\. Equation ( IA3I has the following interpretation: The factor v k i ■ f k i gives the flux of incoming 
particles, while the G- and ^-functions specify the conditions for a collision to take place. A collision between particles k and I 
happens only if the two particles are approaching each other which is ensured by 0(— v k i ■ r k i). At the instant of a collision the 
distance between the two particles has to vanish when two particles touch, which is expressed by S(\r k i | — 2a). Finally, (bf — 1) 
generates the change of linear and angular momenta according to Eqs. (|3} and (|4}. 
The ensemble average, of a dynamic variable, A, is defined by 



(A) t = J dTp(0)A(t) = J dT P (t)A(p) 

n (A4) 
= / \{{d 2 r k d 2 v k duj k ) p(t)A(0) . 



fe=i 



Here pit) = cxp {—iC' + t) p(0) is the iV-particle distribution function, whose time development is governed by the adjoint C* + 
of the time evolution operator £ + . Differentiating equation (IA4> with respect to time yields 

j f {A)t = J dTp{Q)j t A(t) = J dTp(0)iC + A(t) 

= J dTp(0)exp(iC+t)iC + A(0) (A5) 
= J dTp(t)iC+A(0) = {iC+A) t . 



20 



The observables of interest are the averaged energies per particle, or, more specifically, the granular temperatures for the 
two-dimensional system 



fe=i 




(A6) 



k=l 



and the total kinetic energy E = E tr + E rot . To make the temperatures dimensionless we may choose to measure mass in units 
of the particle mass, and velocities in units of the driving velocity vo defined in Eq. 0. 

Assuming a homogeneous density distribution and Gaussian velocity distributions the TV -particle distribution function is given 

by 



""•n, 9( -'-""HiS)^)}' (A7) 

where the product of Heaviside functions accounts for the excluded volume. Hence we get two coupled differential equations 
for the time evolution of the translational and rotational energies 



d , , d 
Jt T ^ =dt 

1 d _ d 

2 <ft ~Jt 



Ttr(t) =Zu(E*)t = (i£+Etr)t = H dr + (i£ + E tr ) t 

T m t(t) =^(E Iot ) t = (i£+E Ioi ) t = (iC + E Iot ) t . (A8) 



The averages on the right hand sides can be calculated as follows. These calculations are almost identical for the translational 
and rotational energies, so we will show in detail the time derivative of 3tr(i) only. 



k^l 

N N . N 



k=l 1=1 " j=l 

iTf(K| 2 + h| 2 ) 

We have used that the binary collision operator iT¥ yields zero acting on any variable other than the ones of the two particles 
involved in the collision. Defining 

N 

dv '■ = n <p r i d2v ] duj j n ©cm - 2a) x 

W (A10) 

/ N N I \ 

eXP ( v -g2^ K|2 -g27-W l ^ |2 J 
and using the definition of iT^ 2 we can write 

N — 1 f 

{i£ + E tI ) t = ~YfdY J dT ' Vl2 ^ ^Cl^is | - 2a) x 

e(-f 12 .* 12 )(6^-i)|(|^| 2 + | W2 | 2 ) 

The change of energy AE tI := y(^+ 2 — 1 ) ( |f 1 1 2 + |f2| 2 ) that results from a collision of particle 1 and 2 depends only on 
the phase space variables of particle 1 and 2. Since we assume spatial homogeneity this change of energy can only depend 
on the relative distance vector r 12 := ri — r 2 as well as the relative translational and rotational velocities v 12 := Vi — v 2 
and u>i2 := u>i + cl> 2 . Further, we assume instantaneous collisions. Therefore the change of energy can only depend on the 
direction of the distance vector fi 2 = • Now we can perform the integrations over those particles that are not involved 
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in the collision. The integrals over d 2 v$ ■ ■ ■ cPv^ and d 2 uj% . . . cPu>n are simple Gaussians. To integrate over d 2 r^ . . . cPriy we 
introduce two more two-dimensional integrals 



J d 2 R 1 d 2 R 2 5 2 (R 1 - ri )5 2 (R 2 - r 2 ) 



over two-dimensional 5 functions, S 2 (r) := S(r x )S(r y ), Using the definition of the pair correlation function 

g(\Ri-R*\) = 
V 2 

N 

I II d 2 r, n ©(l^l - 2a)S 2 (R 1 - ri )S 2 (R 2 - r 2 ) 

j=i m 

N ' 

SU^rjU Q(\ rjl \ -2a) 

where V is the area of the system, we obtain 

iV - 1 / m \ 2 / / \ 
J d 2 R 1 d 2 R 2 d 2 v 1 d 2 v 2 d 2 uj 1 duj 2 g{\R 12 \)x 

exp ^ TUT) ) x 

(H12 • v 12 ) 5{\R l2 \ - 2a)6(-R 12 ■ v l2 ) AE tI . 
Since the change of energy AE tI depends only on R\ 2 := Ri — R 2 , Vi 2 , and u>i 2 , we introduce the variables 

Vi -v 2 U>i+ u> 2 



r := Ri - R 2 , v := 



V2 ' ■ V2 



R := Ri. V .= = — , ft := = — . (All) 

The Jacobian of this transformation is 1. The expression to integrate over is independent of R such that integration over d 2 R 
yields the area V. We write r in polar coordinates (r, </>) and can integrate over dr. Then, we choose the coordinate system for 
integrations over d 2 v such that the unit vector r points in the y-direction. That means we can replace f by the unit vector in the 
y-direction e y and integrate over d<fi which simply yields 2ir. For readability, we use now the unit vector n instead of e y . The 
integrals over d 2 Vand dtt are Gaussians, so that we obtain 



iC + E ir )t = —2-K\p2 a no g(2a) 



I 



2nT tI (t)J \2irT TOt (t) 



(n • v) Q(-n ■ v) AE tr , 

with the number density no := (N — l)/V ?« N/V. 

To solve the integrals above we need to take a look at the change of energy 

-A£ tr := (b\ 2 - l)(\ Vl \ 2 + \v 2 \ 2 ) 
m 

.1 |2 , |„./|2 



= (w^ + \v' 2 n-(\v 1 \ 2 + \v 2 n 

= 477(77 -l)(M 2 -(n- w) 2 ) 

- (1 - r 2 )(n ■ v f + (2an) 2 \h x u\ 2 
+ 4an(2n — 1) (n X oj) ■ v 
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with r) and 7*4(7) given by Eqs. dl9l and 0. Keep in mind that v and uj have been defined as v := Vvi/\2 and u> := U12/V2. 
The difference in calculation for models D and E comes into play at this step [For model D, at this point of the calculation, we 
would express v in polar coordinates (v, 712) and insert 77 = 77(712) instead ofn = 77(7) "assuming" that 7 sa 712. That way 
all integrals become Gaussians and can easily be solved. In particular, the integrals over the last term in Eq. \A12\ vanish.], 
however, we will now go on with model E. To perform the integrations over v and u> we substitute 



w_l ;=v / 2(nxw) = (V2cj,0,0) 
g := \/2 (v + a (n x w)) , 



(A13) 



thus introducing the relative velocity of the contact point g as defined in Eq. 0. The vector lj± points in the x-direction (due to 
our choice of coordinates) and = |u>i2|, so that v xo t = au>±. The Jacobian of this transformation is 2~^. In terms of these 
new variables AE tI reads 



-AE tI = 277(77- l)(\g\ 2 -(h-gf) 
to 



l-r 2 



(n ■ gf + 2ang- u>± , 



(A14) 



and we get 



(iC + E tI ) t 



a n g(2a) 



I 



J d 2 g doj± (h ■ g) 0(— n ■ g) exp 



quia 
4T rot (t) 



(A15) 



exp 



[Isl 



2a £/ ■ u>j_ + a 2 |u)j 



4T tr (t) 



A£ t - 



Next, we express g in polar coordinates (g, 7) where 7 is 770? the usual angle between g and e x but instead - as needed for 
incorporating Coulomb friction - the angle between g and n (i.e. the angle between g and e y , i.e., g = (—gs'mj,gcos'y)). 
Expression (IA14I reads now 



—AE tl = 277(77 - l)g 2 sin 2 7 - ^— - 
to 2 

— 2anguj± sin 7 , 
(note that 7/ = 77 (cot 7) in the Coulomb friction case) and 



-g cos 7 



(A16) 



'iC,E t ,) t = -■ 



(7l- + £/t r /t 

3, 



a n g (2a) 



4T tr (t)/ V4T rot (t) 



OO 



d-l dg / dw^ g cos 7 exp — 



exp 



4T tr (t) 



4T rot (t) 

[g 2 + 2agu± sin 7 + a 2 w 2 ] ) A-Et 



Now we define .A := [TOa 2 /4][l/T tr (?j) + q/T TOt (t)], B := ma sin 7/ [4AT tr (i)] and substitute p 
This leads to Gaussian integrals overp and g. Using x 2 := 1 + T t° t m , we obtain 



2a n g(2a) T t \{t) x A x 



(/-;. 



cos 7 



(1 + (a; 2 — 1) cos 2 7)5 



A (wj_ + So) for wj 



(A17) 



477 



1 

ar 



sin 2 7 — (1 — r 2 ) cos 2 7 
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Up to this point we have not specified, whether we are going to use constant coefficients of restitution or Coulomb friction. All 
this is hidden in r\ which is either a constant or a function of 7. Since we are interested in the Coulomb friction case, we use 
7] = 77(7). We introduce the notation r\ = i#"/imin{| cot 70 1, | cot7|} = minj^o, ^-^l cot7|}, and obtain 



(iC+ E tT ) t = -\J— 2a n g{2a)T^{t) x 4 (H 1 



2 
,.2 



2 V m " u 3 ' ' tr w {3 x 4 
COS 7 



(1 + (x 2 -1) cos 2 7) 5 /2 



1 + r 

2/x — — sin7C0S7 + u 2 (l + r) 2 cos 2 7 



1 \ / cos 7 sin 2 7 

+4?7o 170 7 / d "fl 



(A18) 



x 2 J J '(l + (x 2 -l)cos 2 7) 5 / 2 



70 

After performing the last integration the result can be written in the form 

d 



T tr (f) = H dr -GT?/ 2 {A r (A19) 
at 



770 1 - r] x z 770 x 2 cot 2 70 



2 (l + a; 2 cot 2 7 ) 3 / 2 2 (1 + x 2 cot 2 7o ) 3 / 2 
l+|x 2 cot 2 70 
(1 + x 2 cot 2 70) 3 / 2 



- rj tan 70 1 - 



where G = a uq g(2a), which is the same as Eq. G3> and A r = ±=£-. Similarly T, ot can be calculated using, instead of 

AE tI , the change of rotational energy at collision, 

~AE TOt := (^ 2 -l)(|a;i| 2 + |a; 2 | 2 ) 

= (|a;' 1 | 2 + |^| 2 )-(|a; 1 | 2 + |a; 2 | 2 ) 

m ^^\nxv^ + 4(^-l)(\^-(n.^ (A20) 



a*q" q \q 

aq 



— (2- - 1 I (n x ») • u , 



which can be reformulated as 



- AE mt = ^- (\g\ 2 - (ft ■ g) 2 ) - 2a V g ■ u>± , (A21) 

m q v ' 

using the notation introduced in Eqs. iAl 1> and JA13i . The calculation for the rotational temperature is identical to the one for 
the translational temperature just shown until Eq. JA1 51 into which we insert AE Iot from Eq. 1A21I instead of AE tr . Performing 
the integrals yields 

55 r -®- a ^{(IT?=P^ x (A22) 

^-(x 2 -l)|(l+x J cot 2 7 „) 
?7otan 2 7o / l+|x 2 cot 2 7o 



+ 



{1 + x 2 cot 2 7o) 3 / 2 



Finally, from Eqs. (IA19I > and JA22t the conversant reader may reproduce the transformation to the more convenient coefficients 
in Eqs. ( l33l . 
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For comparison, we quote the equivalent results in three dimensions II 311 : 

3d -T tr (t) = H dl - G 3D T^ (A r - V_o^ML 



2 at i 2 1 + x z cot 70 

r/o /arctan (x cot 70) 1 
2 \ a; cot 70 1 + x 2 cot 2 70 



and 



3 d 



^-77 Trot (0 = G 3D Tj < 770 



3/2 I V 1 



1 _ ( 1 _ 20 \ x 2 



Jot 1 + X z cot 70 

770 , 2 /arctan (a; cot 70) 1 \ 

2 ^ V x cot 70 1 + x 1 cot 2 70 / 



where G 3D = 32 a 2 n .g 3 n(2a) and g 3D {2a) ~ 



